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NOMENCLATURE 


C coefficients in turbulence models 

Cf friction coefficient 

d, D diameter or width of the channel 

E friction factor in law of the wall 

f, f' mean and fluctuating components of variable f 

G buoyancy production/destruction of k = -0 g- u-0 

Gjj buoyancy production/destruction of u]uj = -j3 (gjUj <p + gju^ 0) 

g gravitational acceleration 

H step height 

k turbulence kinetic energy, 1/2 u]uj 

1 mixing length 

1 turbulence length scale 

Nu Nusselt number 

p pressure 

Pe Peclet number 

Ry Reynolds stress 

Re Reynolds number 

R c mean radius of curvature 

R t turbulence Reynolds number 

Rj gradient Richardson number 

Rf flux Richardson number 

Rj c curvature Richardson number 

S swirl number 

U, V, W mean velocity in x, y, z directions 

u, v, w fluctuating velocities in x, y, z directions 

Uj mean velocity component in Xj direction 

Uj fluctuating velocity components in Xj direction 

iv 


L 


U r friction velocity 

x,y,z Cartesian coordinates 

Xj coordinate in tensor notation 

reattachment length 

y + dimensional distance, y \J T /v 

P volumetric expansion coefficient 

T circulation, equation (53) 

Tf false diffusion 

5jj Kronecker delta 

e dissipation rate of turbulence energy 

T) subgrid scale Reynolds stress 

6 polar coordinate 

k von Karman constant 

p dynamic viscosity 

M e ff effective dynamic viscosity (ju + 

v kinematic viscosity 

turbulence kinematic viscosity 
p fluid density 

o^, o e turbulent Prandtl number for diffusion of k and e 
r turbulence shear stress 

r w wall shear stress 

0 fluctuating scalar quantity 

0 stream function 

co vorticity (eqs. 53, 55, and 56), time-mean square vorticity fluctuations (eqs. 7 and 10) 
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TECHNICAL PAPER 


A CRITICAL EVALUATION OF VARIOUS TURBULENCE MODELS 
AS APPLIED TO INTERNAL FLUID FLOWS 

I. INTRODUCTION 


Turbulence is one of the unsolved problems in the area of physical sciences. It is believed that 
the solution of time-dependent three-dimensional Navier-Stokes equations can describe turbulent flows 
completely. However, the computers are not large and fast enough to solve the equations directly, for 
the required range of length and time scales, even for simple flows. Many industrially important flows, 
such as the flow in a space shuttle main engine, are quite complex. Hence, it is of practical importance to 
describe turbulent motion in terms of time averaged quantities rather than instantaneous. This kind of descrip- 
tion leads to the well known turbulence closure problem [ 1 ] . Substitution of apparent mean (Reynolds) 
stresses for the actual transfer of momentum by the velocity fluctuations increases the number of 
unknowns above the number of equations. The problem is then to supply the information missing from 
the time-averaged equations by formulating a model to describe some or all of the six independent 
Reynolds stresses, -p u^Uj. The exact Reynolds stress transport equations can be derived from the time- 

dependent Navier-Stokes equations [2]. These equations express the conservation of each Reynolds stress 
as the Navier-Stokes equations express the conservation of each component of momentum. In turbulence 
modeling one uses a finite number of Reynolds stress transport equations with the supply of missing 
information from experimental results. The time-averaged scalar transport equation contains the turbulent 
heat or mass flux, -p uj0, where 4 > is the fluctuating scalar quantity. One has to model this term to solve 

the scalar transport equation. 

First order closure or mean velocity field closure expresses the Reynolds stresses as a function of 
the mean velocity field. Closure of the Reynolds stress transport equations requires third-order mean 
products of velocity fluctuations to be expressed in terms of second-order mean products (the Reynolds 
stresses themselves) and is known as second order closure or Reynolds stress closure. 

First, this report describes the equations of turbulent flow, and then the turbulence models which 
are often used in engineering calculations. Then it discusses the application of the models and their per- 
formance in internal flows which are of particular interest here. Internal flow encompasses a wide range 
of flows and can be extremely complex. It includes fluid flow in pipes, passages, ducts, conduits, and in 
components such as bends, diffusers, and heat exchangers. Like a majority of all turbulence researches, 
the report deals with incompressible flows. 


II. EQUATIONS OF MOTION 


If the flow variables are assumed to be of the form f = f + f' where f is the mean value of f 
and f' is the fluctuation about the mean, one can obtain the time averaged continuity and momentum 
equations [2,3] 


( 1 ) 



( 2 ) 


)u i _ 9u i _ 1 9p_ m 2 — _ __J 

9t J 3xj p 9xj p ^ * 9xj 


complicated 
the variable 
The various 

turbulence models employed for the closure of the mean equations are discussed next. 


where Ry = u^uj , the Reynolds stress. In order to establish a turbulence model that is less 
than the full equations, it is necessary to establish a way of evaluating Ry in terms of 
u[ , higher-order turbulence statistics and some constants evaluated from experimental data. 


III. TURBULENCE MODELS 


The numerical procedures which are associated with turbulence models to make complete calcula- 
tion methods can be divided into integral and differential types. Differential methods involve direct 
assumptions for the Reynolds stresses at a point and seek the solution of the governing equation in their 
partial differential equation form. 

Integral methods involve the integral parameters of the shear layer momentum thickness, shape 
parameter, skin friction coefficient, etc. One solves a system of ordinary differential equations (for two- 
dimensional flows), whose dependent variables are the profile parameters and independent variable is x; 
in three-dimensional flows, the equations are the partial differential equations in the plane of the layer. 
The important distinction between calculation methods is the type of turbulence model rather than the 
type of numerical procedure. 

The advantage of differential methods is that the restrictions and inaccuracy that arise from the 
need to parameterize the velocity profiles are avoided. Differential methods introduce substantially more 
detailed information about turbulence. Here modeled forms of governing equations and the correspond- 
ing closures for the differential methods are described. 

The turbulence models can be classified in several ways. The one most often used is that arranged 
in order of the number of differential equations solved in addition to the mean flow equations 14] 

(I) Zero equation models 

(II) One equation models 

(III) Two equation models 

(IV) Stress equation models 

Most of the models, classes (I) to (III), use Boussinesq eddy viscosity model. Bradshaw, et al. [5], 
however, assume the constancy of the r/pk ratio, in boundary layer flows. Here r is the shear stress and 
k is the turbulent kinetic energy. It is important to note here that in this case the mean momentum and 
continuity equations form a hyperbolic system in contrast to the parabolic system obtained with the use 
of eddy viscosity models [6], Other models which do not use the eddy viscosity assumption (class IV) 
obtain the Reynolds stress from a differential equation. 


Zero equation model, which uses only the partial differential equation for the mean flow field 
and no transport equations for turbulence quantities, is also called “mean field” closure [7], The classes 
(II) to (IV) are called “transport equation” closures. Yet another classification is based on the highest 
order of velocity product for which a transport equation is used. Zero equation models use partial differ- 
ential equations for the Uj only and are therefore “first order” models. Classes (II) to (IV) use partial 

differential equations for ujuj and are “second order” closures, while some class (IV) models use transport 
equation for third order products u-UjUj [4]. 

Bradshaw [8] describes the interplay between the development of models and the experiments. 
The monograph by Launder and Spalding [9] gives the mathematical concepts of turbulence models. 
Bradshaw and Cebeci [10] describe the calculation methods for various classes of turbulent flows. 
Gosman, et al. [11] present various aspects of computation of recirculating flows. They form the basis 
of many of the current turbulent flow computational schemes. Several reviews have appeared recently 
concentrating on specific aspects of turbulence modeling [12 to 21]. This report concentrates on the 
turbulence models and their applications to internal flows. In the following, the main classes of turbu- 
lence models are described briefly. 


3.1. Zero Equation Models 

Zero equation models are mostly based on the eddy viscosity concept. This concept stems from 
the convenience associated with maintaining an approach, for turbulent flows, which is similar to that 
for laminar flows. This convenience of concept is coupled with the possible mathematical convenience 
of retaining the same form of differential equations for laminar and turbulent flows and allowing the 
use of the same solution procedure. 

The first turbulence model proposed, the Prandtl’s mixing length hypothesis, is still among the 
most widely used models. It employs the eddy viscosity concept which relates the turbulent transport 
terms to the local gradient of mean flow quantities. For example, for thin shear layers 


- uv = 



(3) 


where 


i> t = eddy viscosity. 


The Prandtl mixing length hypothesis calculates the distribution of eddy viscosity by relating it to the 
mean velocity gradient 


*1 = c p V 


9U 

9y 


(4) 
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This relation involves a single unknown parameter, the mixing length l m whose distribution 
over the flow field has to be prescribed with the aid of empirical information. is a constant. The 

mixing length model has been used for thin shear layers and wall boundary layers [16]. The main draw- 
back with this model is the evaluation of l m for different flows. The evaluation of l m becomes difficult 

for recirculating flows and three-dimensional flows. In the already empirical specification of the mixing 
length, it is difficult to incorporate in any useful manner, the effects of curvature, buoyancy or rotation. 
The transport and history effects of turbulence are not accounted in the mixing length model. The more 
generally applicable models to be described below account for these effects by introducing transport 
equations for turbulent quantities. 

3.2. One Equation Models 

The one equation model requires the solution of an equation for the turbulent kinetic energy, k, 
and, as a result, allows for its transport. The turbulent kinetic energy equation can be derived from the 

Navier-Stokes equations [12]. The eddy viscosity is modeled by k^ 2 1. The length scale 1 is 

specified algebraically and hence is flow dependent. This approach is not very popular since it performs 
only marginally better than the zero equation model. 


3.3. Two Equation Models 

This class of models is the one widely used in present day engineering calculations. In attempts 
to eliminate the need for specifying the turbulence length scale as a function of position throughout the 
flow, a second differential equation which in effect gives 1 has been used. 

In general, one looks for an equation for a quantity that is a combination of k and 1, Z = k a 1^. 
Such an equation can also be derived from the Navier-Stokes equations. It has the form [12,22], 


P 


DZ 

dT 


9 / ^t 9z\ 

9y \a z 9y J 


+ Z 




(5) 


Here a z is a Prandtl number for the diffusion ofZ,S z is a secondary source term which appears in some 
models, and Cj and C 2 are constants. 

The Imperial College group lead by Prof. Spalding has experimented with three different kinds 
of two equation models [21]: k-kl; k-co, and k-e. 1 is a length representing the macroscale of tur- 

bulence which may be defined in terms of k, e and a constant Cq through 


1 = C D k 3/2 /e 


( 6 ) 


'l 

co is a quantity having the dimensions of (time)" - which has been interpreted as representing the time 
average square of the vorticity fluctuations, and it can also be defined in terms of k, e and Cq. 
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(7) 


co = e 2 /(C D k) 2 


and e is defined by 


3U| 0U| 

3x k 9* k 


( 8 ) 


The above definitions imply that 


dkl dk de 

kT ~ 2 T ~ 7" 


dto 

co 



(9) 

( 10 ) 


With the aid of these equations, it is possible to transform one pair of equations into another [21]. Thus, 
the three models differ only in the mathematical form but not in content. However, the k-e model has 
become most popular because of the practical advantage that the e-equation requires no extra terms near 
walls. Also, e itself appears in the k-equation and the e equation requires no secondary source term. 
Hence, only the k-e model is described in detail here. 


The k-e model employs the eddy viscosity and relates the eddy viscosity to k and e. One solves 
two differential equations, one for turbulent kinetic energy and the other for its dissipation. The modeled 
equations for k and e are given below [21 ] : 


(a) Kinetic energy equation: 


Dk _ 1 9 "Mt 9k + Mt / 9u i + 9U k\ 9U i 

Dt p 9x k o k 3x k p \ 9x k 9xj J 9x k 


(b) Kinetic energy dissipation rate equation: 


De 

Dt 


I JL 

P 9x k 


^t 3e 
~ 9x k _ 


C 1 H e 
P k 


\9x k + 9 Xi 



9x k 



( 12 ) 


Mt ~ Cju P k 2 /e 


(13) 
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The constants in these equations have been found to take the following values [21]. 


C = 0.09 
Cj = 1.44 
C 2 = 1.92 
a k =1.0 
a e = 1.3. 


The coefficients are constants in the sense that they are not changed in any calculation. However, 
these constants need to be changed in order to accommodate the effects such as curvature, low Reynolds 
number, near wall, etc. 

Close to the walls, the local Reynolds number of turbulence is small and the molecular transport 
becomes important. To resolve the wall layer properly in a numerical solution procedure one requires an 
extremely fine mesh near the wall. Fine mesh size increases the computational time by an order of 
magnitude. Hence, most computational schemes avoid it by using the wall function to bridge the whole 
of the wall layer to the fully turbulent region. The first computational point p near the wall (Fig. 1 ) is 
to be located in the fully turbulent log-law region. The momentum flux at such a point p satisfies the 
relation [21] 


U. 


— Cj/ 4 k 1/2 = I i„ | E y 
P K 




(T/P) w 


(14) 


Here Up is the velocity of the fluid at the point p along the wall and r w is the shear stress at the wall 
in the direction of Up. The quantity kp, the value of k at the grid point P, is computed by assuming 

that the generation and dissipation of energy are equal in the wall layer where the shear stress is uniform 
and the length scale is proportional to the distance from the wall. The quantity E is a function of wall 
roughness and is equal to 9 for smooth walls. Extensive accounts of wall functions and their manner of 
application are discussed in References 11 and 21. Recent formulations of near-wall treatments are 
discussed in Section IV. 
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3.4 Reynolds Stress Models 

The two equation models which have been extensively used in engineering calculations have 
several limitations. The main practical limitation is the assumption of isotropic eddy viscosity. The same 
values of i> t are taken for different' u^uj terms. The effects of curvature, rotation and buoyancy forces 

have to be modeled separately. In the models discussed so far the local state of turbulence is assumed 
to be characterized by one velocity scale -^/k. The individual u-Uj is related to this scale. In actual flows, 

the scales may develop quite differently. In order to account for the different development of the 
individual stresses, transport equations for u|uj have been introduced. These equations can be derived in 

exact forms but they contain higher order correlations that have to be approximated in order to obtain 
a closed system. In this model one needs to solve the equation for the turbulence energy dissipation rate 
e, in addition to those for ujiH for the length scale. A particular advantage of the Reynolds stress model 

is that terms accounting for buoyancy, rotation and other effects are in principle introduced 
automatically. However, many problems arise in the solutions of the model equations as discussed below 
and in Section 3.4.1. 

Models employing transport equations for u-uj are called second order closure models. Several 
closure schemes have been proposed. The well tested one is that of Launder, et al. [23] . 

Launder et al. proposed the following model transport equations for u^uj 



Convective Diffusive Pjj = stress production 

Transport Transport 



Pressure Strain 


- (3 (gjUj </> + gjUj 0) - j e 5y (15) 

Viscous Dissipation 

Gjj = buoyancy production. 


The rate of change, convective transport, and mean field as well as buoyancy production terms 
are exact whereas the diffusion, pressure-strain/scrambling, and viscous dissipation terms are model 
approximations. The diffusion fluxes of UjUj have been expressed by simple gradient diffusion models. 

Local isotropy has been assumed to prevail so that the dissipation is the same for all the three normal 
components (and thus 2/3 of the total dissipation e) and, so that the viscous destruction terms for the 
shear stresses are zero. The most important assumption concerns the pressure-strain/scram bling terms, 
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since for shear stresses these are the main terms to balance the production of these quantities. The 
pressure-strain/scrambling model consists of three parts. The first one represents the interaction of 
fluctuating components only, the second the interaction of mean strain and fluctuating quantities and the 
third the effect of buoyancy forces. Several versions of the pressure strain model have been proposed to 
correctly predict the experimentally observed results. To account for the wall damping effects a wall 
correction must be introduced in the pressure-strain model. Launder et al. make the empirical constants 

in the pressure-strain model a function of the relative distance from the wall, 1/y a k^/fey). Because 
of the complexity and the large amount of computational efforts required, the model has not been as 
widely used as one would like it to be. 

3.4.1 Invariance and Realizability of Reynolds Stress Models 

Turbulence models should be broad based so that they can be applied to most practical problems 
without “tuning.” The Reynolds stress models come close to this ideal, since they account for more 
physical processes such as those due to curvature, rotation, buoyancy, etc. automatically. For general 
applicability, these models have to satisfy many constraints [17]* Of these constraints, tensor invariance 
and realizability seem to be important. The first one requires that the replaced terms have the same tensor 
form as the original terms. This will ensure that they transform properly in different coordinate systems. 
Such type of modeling is called invariant modeling [20]. The second constraint, realizability was first 
introduced by Schumann [24]. This condition requires that the equation for turbulent stresses have 
the property that all component energies remain non-negative and all off-diagonal components of 
Reynolds stress tensor satisfy Schwartz’s inequality. These and the additional condition put forth by 
Schumann are written in a numerically most convenient form [24] : 


RH>0 


(16) 


R 11 R 22" R 1 2 2 >0 ’ 07) 

R 1 1 (R 22 R 33 " R 23 ) ‘ R 12 (R 12 R 33 " R 23 R 13 ) + R 13 (R 12 R 23 “ R 22 R 13 ) > 0 • (18) 


Similar conditions apply for scalar fluxes. Schumann shows that the exact Reynolds stress 
equations satisfy the realizability condition. He also indicates the nonrealizability of some of the existing 
models. However, attempts to satisfy realizability conditions lead to invariably complicated model 
expressions [14]. Hence, at present no special efforts are taken to see that model equations satisfy 
realizability conditions. Care must be exercised in implementing and interpreting these models. 


3.5 Algebraic Stress Model 

In Reynolds stress models, there are differential transport equations for each component of ujUj 

in addition to e equation. To reduce the computational effort, Rodi [25] proposed an algebraic relation 
for calculating the Reynolds stresses. The convection and diffusion terms in the transport equations of 
UjUj are replaced by model approximations, reducing the equations to algebraic equations. Rodi assumes 

that the transport of uJTTj is proportional to the transport of k and that the proportionality factor is 
UjUj/k. With this approximation incorporated, the transport equations yield algebraic expressions for 
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UjUj that contain the various production terms appearing in the u-m equations. Thus the gradient of 

mean flow quantities, k and e appear also in the expressions, so that k and e equations have to be added 
in order to complete the turbulence model. The algebraic expressions together with k and e equations 
form an extended k-e model. 

Algebraic stress models are suitable whenever the transport of ujuj is not important. Algebraic 

stress relations are basically like eddy viscosity formulations [15] and therefore are not applicable to 
cases where countergradient transports occur (Section VII). On the other hand, all effects that enter the 
transport equations for ujuj through the source terms for example, body force effects (buoyancy, rotation 

and streamline curvature), nonisotropic strain fields and wall damping influence can be incorporated. 
Algebraic stress models can therefore also simulate many of the flow phenomena that were described 
successfully by stress equation models. 


3.6 Multiple-Scale Models 

All the turbulence models discussed so far are based on the assumption that in all the flow situa- 
tions each variable has a spectrum of universal form which can be characterized by the scales of the 
energy containing range. Difficulties arise when the spectrum is not an equilibrium one or when the flow 
exhibits distinctly different ranges of scales. Complex models have been devised to predict such flows 
[17]. The model of Hanjalic, et al. [26] appears to be a tractable and useful one. They split the spectrum 
into a large scale part and a small scale part with different time scales for transfer into the large scale 
part and transfer of the large scale part to the small scale part. A brief description of their model is given 
below. 


The turbulence spectrum consists of independent production, inertial and dissipation ranges. 
Kj denotes the wave number above which a significant mean strain production occurs while K .2 is the 

largest wave number at which viscous dissipation of turbulence is unimportant (Fig. 2). Energy leaves the 
first region (production) at a rate e n and enters the high wave number or dissipation region at a rate e. 



Figure 2. The spectral division for multiple scale model [26]. 
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Between the two regions, occupying the intermediate range of wave numbers is the transfer 
region, across which a representative spectral energy transfer rate e-p is assumed. This simplified energy 
spectrum is the basis of Hanjalic et al.’s model. 

The total turbulence energy k is assumed to be divided between the production range kp and the 
transfer range kp. At high Reynolds numbers there is negligible kinetic energy in the dissipation range. 

In a homogeneous flow the levels of kp and ky are controlled by the transport equations, 


Dk p _ aui 

"5T = ' U i U j ~ e P 


(19) 


P 


k 


Dk T 

~dT = € p 


- e 


( 20 ) 


where denotes the production rate of turbulence energy by mean shear which is assumed to be 
contained in the wave numbers below Kj. The following equations have been proposed for characterizing 
the evolution of these transfer rates. 


Dej _ 

Dt Cel 


p kT" 


(C e2 + C e3 A) + D el 


( 21 ) 


The quantity “A” refers to the turbulence anisotropy defined as 


(ujUj - 2/3 Sy k) (u^. - 2/3 5jj k) / k 2 . 


with C e j, C e -> and C e3 are constants determined from experiments. D f j denotes the diffusive transport 
[equation (21) is for single scale model]. Equation (21) makes the local rate of ej dependent on the local 

mean strain rate and the anisotropy of the stress field neither of which under conditions of local isotropy 
can directly affect the dissipation rate. Equation (21) may be regarded as a spectral energy transfer asso- 
ciated with large scale interactions. The proposed equation is: 


D e e p e p 2 

— — — Cpi p l- Cp^ + D,- n 

Dt P1 k k p P2 k p ep 


( 22 ) 


where the partioned energy k p replaces the total energy giving, as the characteristic time scale, the turn- 
over time of the large scale motions. C p j and C p2 are constants. 
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The proposed form of the ey transport equation is 


~dT 




e P e T 

k T 



( 23 ) 


Cy | and Cj2 31-6 constants determined from experiments. Energy dissipation rate now responds only 

slowly to the applied mean strain. This feature makes the form (23) more capable of representing the 
rapidly changing fields than equation (21). Hanjalic et al. computed the self-preserving plane and 
axisymmetric jets with this model using simple turbulent viscosity hypothesis for the shear stress uv. 
They obtained very good agreement of the spread rates of the jets with the experimental values. 

We have given here a brief account of various turbulence models along with their main features. 
Second order models and the models for two phase flows are not considered here. In the next section 
we discuss the corrections made to the widely used two-equation model to extend its applicability. 


IV. MODIFICATIONS TO k-e MODEL 


Of all the models discussed in Section III, the two-equation, k-e model has been widely used to 
predict various classes of flows of engineering importance. It has been used for predicting flows with 
heat/mass transfer, streamline curvature, etc. These effects have been accounted for by correcting the 
model equations. The low Reynolds number modification, wall-layer models and the streamline curvature 
corrections are discussed here. 


4.1 Low-Reynolds Number Model 

Jones and Launder [27] extended the k-e model to low-Reynolds numbers so that the turbulence 
model equations can be valid throughout the laminar, transition and fully turbulent regions. In this 
version of the model, k and e are determined from the following pair of equations. 



In these equations, Cj, and o £ retain the values assigned to them for high Reynolds numbers while 
and C 2 vary with turbulence Reynolds number. 

exp [-2. 5/(1 + R t /50)] , (26) 
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(27) 


C 2 = C 2oo [1-0.3 exp (-R t 2 )] , 


subscript refers to fully turbulent values. We note here, that the laminar diffusive transport becomes 
of increasing importance as the wall is approached and the extra destruction terms included are of some 

significance in the viscous and transitional regions. The term 2.0 i>M t /P (S^Uj/dXjBxj) 2 in the e equation 

produces satisfactory variation of k with distance from the wall. The measurements indicate that the level 
of energy dissipation rate is constant in the immediate neighborhood of the wall. In the computations e 

is set to zero at the wall and an extra term, -2v (3k 1 /^/3 X j) 2 is introduced to the k equation. This extra 
term is exactly equal to the energy dissipation rate in the neighborhood of the wall. 


4.2 Near-Wall Models 

4.2.1 Two-Layer Model of Chieng and Launder 

Detailed modeling of the near wall region is necessary for better prediction of the flow field. 
Recently Chieng and Launder [28] introduced detailed near wall modeling in their computation of flow 
in a pipe expansion. Their model is briefly discussed here. Figure 3 shows a scalar node P whose 
associated volume is bounded on the south side by a wall. The near wall flow is treated as viscous (but 
not laminar) out to a distance y v from the wall and fully turbulent beyond this. The sublayer thickness 

y v is such that the Reynolds number y y k y ^~/ p (= R y ) is constant, taken equal to 20. The near wall 

cell is large enough that the node P lies outside the viscous sublayer. Over the fully turbulent region 
the mean velocity parallel to the wall is assumed to vary with the distance from the wall according to 


i V 


U k '/ 2 /(r w /p) = — In E* 


(28) 



Figure 3. Near-wall two-layer model. 
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An expression for r w 


can be written in terms of U • 

P 


r 


w 


K* P Up ky 1 / 2 

In E* y p ky 1 / 2 /^ 


( 29 ) 


The equations (28) and (29) are the same as the low Re version of Launder and Spalding [21] except 
that the kinetic energy term is evaluated here, at the edge of the viscous sublayer, rather than at node P. 
The k* in equation (28) is proportional to the von Karman constant k. In local equilibrium, where 

- 1/2 

k = (r /p) C aoo the velocity relation should reduce to the conventional logarithmic law: 

W r 



y V( t w /p) 

V 


(30) 


The limiting form is achieved by taking k* = k (0.23) and E = E* The constant E* 

can be evaluated in terms of R y by noting that within the viscous sublayer the velocity displays a linear 
variation which may be written as 


Ukyl/ 2 _ykyl / 2 
(r w /p) v 


(31) 


At y y we find that E* - exp (k* Ry)/R y . Thus for the assumed values of k* and R y , E* = 5.0. 

The linear variation between the node P and its neighbor (point N) is used to extrapolate to the 
viscous sublayer (i.e., k y in terms of kp and kj^). Within the viscous sublayer a parabolic variation of 
kinetic energy is assumed. 


k = k v <y/y v )2 


(32) 


This corresponds to the linear increase of the fluctuating velocity with distance from the wall. 

The dominant contribution to the energy generation rate is the term r-p 9u/9y, the rj is the local 
shear stress. The mean generation rate per unit volume is 


1/y n / 


9U , 

l T w + ^ r n " r w) y / y n-* y 


(33) 
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Using equation (30) for the mean velocity distribution 


T W ( U n -U v) r w ( T n - r w) / 
Mean generation rate = + - I 1 

y n P «.* ky 1 / 2 y n \ 

Unlike the generation, the dissipation rate of turbulence energy 
equal to 2v (9k^ 2 /9y) 2 . Using equation (25) we obtain 


2v k v 



The mean value of the dissipation rate can be obtained by integrating (35) between y y and 



is not zero in the sublayer. 


( 34 ) 

It is constant 


(35) 


4.2.2 Two- and Three-Layer Models of Amano 

Amano's near wall model [29] is the same as that of Chieng and Launder, except for the treat- 
ment of the generation and destruction terms in the e equation. In Reference 28 the value of e in the 
near wall cell was approximated under local equilibrium condition. Amano develops the treatment of 
the e equation near the wall cell, taking into account that the value of e near the wall is an order of 
magnitude larger than that in the fully turbulent core, and reaches its maximum at the wall. Each term 
in the e equation should be evaluated in accordance with the k equation rather than approximated 
under local equilibrium conditions. He has considered both a two layer model and a three layer model. 
In the two layer model the region is divided into two distinct parts: a viscous sublayer region (0 < y + < 
1 1 ) and an overlap region (11 < y + < 400). The e equations for these regions have been developed [29] . 

In the three layer model he tries to approximate the experimental velocity profde more closely 
by dividing the near wall region into three layers: a viscous sublayer (0 < y + < 5) adjacent to the wall; 
a buffer layer (5 < y + < 30); and an overlap layer (30 < y + < 400) (Fig. 4). 



(a) WALL ADJACENT CELL (b) TURBULENT KINETIC (c) ENERGY DISSIPATION (d) TURBULENT SHEAR STRESS 

ENERGY RATE 


Figure 4. Near-wall three-layer model [29], 
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He assumes the following variation of k, e and r in the three regions: 
Viscous Sublayer 0 < y + < 5 


k = k y (y/y v ) 2 

( 9k 1/2 2 
e 2 \ ay 


r = 0 


(36) 


Buffer Layer (5 < y + < 30): 


k = k B y/y B 
e = k 3 / 2 /q y 
r = r B (y/y B ) 3 . 


Fully Turbulent Region (30 < y + < 400): 


, k n - k B / K P 

k = y + Ikp - 

y n - yB V yp 


/ k P -k N \ 

\ p yp - yN yp J 


by + a 


e = k 3 / 2 /C, y 


T = T w + ( T n - r w ) y/y n 


(37) 


(38) 


where 


a = 


kp - 


k P " k N 

y p 

yp - yN 



y n - yB 


He obtains the mean generation and destruction rates for k and e equations, incorporating the above 
relations. He claims to obtain better results with these near wall models, the three layer model being 
superior to the two layer one. 
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The near-wall models of Amano and Chieng and Launder have been useful for improving the 
prediction of wall heat transfer rates. However, it appears that these models do not improve significantly 
the flow field predictions. For example, the models have little effect on the prediction of the location 
of the reattachment point [28] (also, see Section V). 


4.2.3 PSL Approach 

In an attempt to model the physics of near wall region, Iacovides and Launder [30] have recently 
proposed a numerical scheme which employs the classical boundary layer assumption for the thin wall 
layer and eliminates the use of wall functions. In this approach a thin parabolic sublayer (PSL) near the 
wall is assumed across which static pressure variation is either negligible, or if the surface is highly curved, 
can be obtained assuming radial equilibrium. The PSL is taken to extend over the whole low Re region 
where the transport coefficients vary by orders of magnitude. A fine computational grid treatment for 
this layer eliminates the use of wall functions. In PSL approach no pressure need be computed or stored 
at the grid points within the layer and the velocity component normal to the wall at these node points 
are obtained by cell continuity rather than solving the normal momentum equation. Iacovides and 
Launder employ a staggered arrangement of dependent variables as shown in Figure 5, and the normal 
component of velocity evaluated from 


V(I,J) = [U(i-1 ,J) - U(I,J)] 5y/5x + V(I,J-1 ) 


(39) 




PRESSURE NODE 


V - VELOCITY OBTAINED 
FROM Y- MOMENTUM 
EQUATION 


V - VELOCITY OBTAINED 
FROM CONTINUITY 


U - VELOCITY 
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The PSL approach has been successfully used for computations of flow and heat transfer in a 
pipe expansion, in a spirally fluted tube and in a 90 deg tube bend by Launder [31]. It has been 
found that the PSL approach is not only as economical as the wall function approach to obtain grid 
independent solutions, but also produces more realistic solutions. For example, in the flow in a pipe 
expansion, the PSL approach reveals a secondary recirculation region, not produced by the wall function 
method, for high expansion ratios [31]. The secondary recirculation region accounts for the Nusselt 
number variations downstream of the expansion observed in the experiments. 


4.3 Curvature Effects 

An exhaustive review of the subject and many insights into the effects of streamline curvature 
on turbulence and ways of incorporating into the turbulence models were provided by Bradshaw [32]. 
Large changes in Reynolds stress observed in the curved flow experiments cannot be accounted for by 
the extra terms appearing in the mean flow and Reynolds stress transport equations. Streamline curva- 
ture causes large changes in higher order quantities (second order) of the turbulence structure. Effects 
are similar when the streamlines have a component of curvature in the plane of mean shear irrespective 
of the cause of the curvature, caused by surface curvature or swirl or rotation of the whole system. 
Since the models can recognize only the explicit extra terms due to curvature in the mean motion and 
turbulent transport equations, they will severely underpredict the effects. Hence, the models have to be 
modified to account for the curvature effects. 

Modifications to the k-e model have been made to incorporate the streamline curvature caused 
by surface curvature as well as by swirl or rotation. The modification of the length scale suggested by 
Bradshaw is employed in most computations. It is made a linear function of the curvature Richardson 
number. Another way of incorporating the curvature effect is through the modification of C^. For the 

convenience of discussions on model predictions in Section V, the modifications employed in the flows 
with surface curvature and in the flows with swirl are given here separately. 


4.3.1 Flows with Surface Curvature 

Gibson [33] derived a set of algebraic stress equations from which the effect of curvature on 
mixing length and turbulent Prandtl number can be obtained. His analysis was, however, limited to 
moderate curvatures. Following the idea of Gibson, Pourahmadi and Humphrey [34] studied the 
effect of curvature for any arbitrary curvature. They obtain a general expression for by equating the 

algebraic expression derived for ugu r with the Boussinesq approximation for ugu r in which is given 
by 


Ht/P = (k 3 / 2 /e) k 1 / 2 . 


They use the expression derived for in the k-e model equations. The variation of for a curved 
channel obtained by them is shown in Figure 6. 

The figure shows the variations of as a function of radial position in channels of different 
curvature. In general, is seen to increase at both walls of a curved channel at a rate inversely propor- 
tional to the channel curvature (R c /D). For strong curvatures they obtained unrealistic values of the 
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parameter in the region 0.3 < r < 0.65 due to the lack of consistency of the ratio u^Uj/k. They used the 
fixed values of = 0.09 in the region 0.3 < r < 0.65. 




Figure 6. Surface curvature: geometry and variation of for fully developed 
curved and straight channel flow [34]. 

The relative importance of wall curvature and wall pressure fluctuations on CL have been shown 
by two sets of profiles (R c /D = 10 and 20) with imposed (m=0) symmetric wall functions of f 

equivalent to the straight channel flows in so far as the wall pressure corrections are concerned, while 
retaining the direct influence of wall curvature on CL. The C profiles show that, the curvature at the 

p* r 

concave wall acts to enhance while the curvature at the convex wall acts to suppress it. The inclusion 
of wall pressure corrections in the pressure strain (m = 2.56 and m = 7.95) further increases C„ at both 

the walls, but at the convex wall the direct influence of curvature effects ultimately dominate the wall 
pressure contributions to causing a net decrease in its value with increasing distance from the wall. 

Launder et al [35] introduce corrections to the length scale determining e-equation to account 
for curvature effects. The constant C 2 is made to depend on the gradient Richardson number defined 
as 



W 3(rW) 


(40) 


They solve the following form of k-e equations in addition to the mean flow equations for two dimen- 
sional curved duct flows. 
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( 41 ) 



(42) 


where 




The new constant introduced C c for curvature varied from 0 to 0.5. 

The influence of curvature can also be accounted for by modifying the velocity scale. Since 
velocity scale effects are included in k-equation such modifications have not been very useful [36-38]. 
Lilley [39] includes the effects of curvature in kl equation by an additional term in his k-kl two equa- 
tion model. 


4.3.2 Flows With Swirl 

For flows with swirl, the k-e model is varied in a simple way by introducing corrections to the 
length scale determining e-equation. There are two corrections in use. One is to modify the constant 
C 2 appearing in the sink term and the other is to modify the constant Cj appearing in the source term 

of the e-equation, by relating them to the Richardson number. However, the form of the Richardson 
number used in the two cases are different. The term C 2 is corrected using the gradient Richardson 

number Rj given by equation (40) [40,41,35], For example, Launder ei ai. [29] replace C 2 by 


C 2 = 1.92 (1 - 0.2 Rj) 


(43) 


Srinivasan and Mongia[41] split R^ into two parts: (i) swirl Richardson number, equation (40), and 
(ii) curvature Richardson number defined below, in an arbitrary way. They represent C 2 as 


C 2 = 1.92 exp (2^ R t + 2« c R^) 


(44) 


R ic = y* + v 2 /R c 


1 a 

- r- (rU) + 
r 9r 



(45) 
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where the radius of curvature is given by 


(1/R C ) = 


UV [1/r 3/dr (rU) - (9U/9x)] 
(U 2 + V 2 ) 3 / 2 


( 46 ) 


The values of a c and were arrived at by parametric studies, to be a c = -0.75 and = -2.0. 

S L- S L/ 

Rodi [42] is probably the only one employing the correction to Cj by relating it to the flux 
Richardson number Rf 



He represents Cj as 

Cj = 1.44 (1 +0.9 R f ) . (48) 


This form of correction has been found to produce the destabilizing effect of rotation on turbulence 
in swirl flows correctly. The use of Rj- is to be preferred for flows with swirl for the following reasons: 

(i) It represents correctly the effect of rotation on turbulence. 

(ii) It is the only suitable parameter for three dimensional flow analysis. 

(iii) It reduces to gradient Richardson number when the curvature is small. 

The curvature effects may, however, be taken care of automatically, if one uses transport equations for 
turbulent shear stresses. 

The applications of the various turbulence models and the model varients discussed here to 
internal flow predictions are discussed in detail in the following section. 


V. APPLICATIONS OF TURBULENCE MODELS TO THE PREDICTION 

OF INTERNAL FLOWS 


The geometrical configuration, the fluid and the flow rate all determine the complexity of the 
flow. For any combination ot these factors, one is interested in a prediction method to obtain informa- 
tion on pressure drop, flow pattern, temperature distribution, etc. This paper is concerned here with 
the computation of turbulent flows via the solution of time averaged Navier- Stokes equations tor a set 
of flows ol engineering interest. Turbulence closure models and their prediction capabilities will be 
considered. In the 1 981-Stantord conference on complex turbulent flows [43] several isothermal flows 
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of engineering interest were computed and compared with carefully compiled experimental data. This 
included a number of confined flows, which they classified into two groups: (i) attached flows and 
(ii) separated flows. Here, the application of turbulent closure models to a selected set of flows of 
interest is discussed and the performance of the models is evaluated. 

Precise comparison of the performance of closure models is difficult because of the difference 
in numerical schemes used. Turbulence models and numerical methods are considered independent ele- 
ments of the total turbulent flow computation. In principle, one can solve any set of turbulence model 
equations with any numerical method. In practice, when the agreement between the computed and 
experimental results is not so good, one is not sure where to put the blame; i.e., turbulence model or the 
numerical method; many tend to blame the turbulence model, because of the faith in their numerical 
method. Hence, an evaluation of the turbulence closure models rests mainly on their predicted flow 
quantities of engineering interest (such as the size of the recirculating flow, profiles of velocity, pressure, 
kinetic energy, etc.) and their closeness to the experimental or expected results. 

In this section, an attempt is made to catalog recent applications of turbulence models to con- 
fined flows under the following classes: 

5.1 Plane two-dimensional computations 

(a) Backward facing step flow 

(b) Symmetric expansion 

(c) Flow over a square obstacle. 

5.2 Axisymmetric two-dimensional computations 

(a) Flow in a sudden pipe expansion 

(b) Diffuser flows 

(c) Flows with curvature and swirl 

5.3 Three-dimensional computations 

Non-circular duct flows 

In each class, the predictions by various turbulence models are compiled and the major predicted 
quantities are compared. Some observed discrepancies between the prediction and experiment are 
discussed. 


5.1 Plane Two-Dimensional Flows 

5.1.1 Backward Facing Step Flow 

The flow over a backward facing step (Fig. 7) was a bench mark problem in the 1981-Stanford 
Conference: flow No. 0421. The results of different turbulent models [44-50] are briefly reviewed here. 
Since we are interested in the recirculating flow, some special methods (which do not solve an elliptic 
equation) adopted for this problem are not included here. Recent results of Ilegbusi and Spalding [51], 
Nallasamy [52] and Syed et al. [53] are also included. Table 1 lists the turbulence models used by 
different authors and some observations. 


21 





Figure 7. Flow over a backward facing step: predicted separation length, X^. 

A useful parameter for comparison in backward facing step flow is the length of the recirculation 
region x^. Figure 7 gives the reattachment lengths predicted by different computations. The experi- 
mentally observed value of x^ is (7.0±0.5) step height [54]. Standard k-e models underpredict the 

reattachment length by as much as 20 percent. Modified k-e models [49,51], algebraic stress model [47] 
and Reynolds stress model [50] predict recirculation length fairly well. Demirdzic et al. [49] use the 
near wall model (Section 4.2) of Chieng and Launder [35] with the standard k-e model, with modifica- 
tions for the production and diffusion terms in the e equation. This computation underpredicts the value 
of xr by about 1 1 percent. It is observed that the better prediction of x^ is due to the modification of 

the diffusion and production terms to account for the streamline curvature. The recent study of Ilegbusi 
and Spalding [51], using the near wall treatment [28] and modified k-w model, predicts x^ which agrees 

with the experiments [54,55], This model requires the specification of two additional constants. The 
exact reason for the underprediction of x^ by the standard model, however, is unclear (see Section 
5.2.1). 


All the methods overpredict the shear stress in the separated shear layer. This implies that in this 
region, the computed turbulent viscosity will be higher than that existing in an actual flow. This, of 
course, can result in a shorter reattachment length as observed in all the predictions. The redevelopment 
of the flow downstream of the reattachment is also not correctly predicted in any ot the predictions. 
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TABLE 1. FLOW OVER A BACKWARD FACING STEP 


Reference 

Turbulence Model 

Flow Properties 
Compared 

Remarks 

Mansour and Morel [44] 

k-e 

Wall static pressure P 
Mean velocity U 
Shear stress uv 

x r =5.2H 

Different near wall models 
have no effect on x R 

Pollard [45] 

k-e 

p, U, uv 

x R — 5.88 H 

Reasonable agreement with 
cross stream uv, poor predic- 
tion of axial variation of uv 

Rodi et al. [46] 

k-e 

p, U, uv 

x R = 5.8 H 

Shear stress over predicted, 
recovery beyond reattachment 
not well predicted. 

Launder et al. [47] 

Algebraic stress model 

p, U, uv 

X R “ 6-9 H 

Agrees with the experimental 
result. Predicts a slower down- 
stream development compared 
to experiment. 

Spalding et al. [48] 

k-e 

p, U, uv 

x r = 6H 

Demirdzic et al. [49] 

Modified k-e. Modifica- 
tions to production and 
diffusion terms in e 
equation. Near wall model 
of C and L [28] 

p, U, uv 

x R = 6.2 H 

rate of growth of outer edge 
underpredicted. Effect of 
each of the two modifica- 
tions is not known. 

Donaldson et al. [50] 

Reynolds stress model 

p, U, uv 

x R = 6. 1 H 

Prediction is good but 
expensive. 

Ilegbusi and Spalding [51 ] 

k-e modified 
near wall model of 
C and L 
k-co modified, 
new term added (two 
more constants) near 
wall model of C and L 

p, U, uv 

x R = 7.2H 

k-e and k-w predict equally 
well. Predicts the shear stress 
variation across the channel 
as well as along the channel 
reasonably accurately. 

Nallasamy [52] 

k-e model 

P> U, uv 

Contours of k, e 

X 

n 

oo 

Syed et al. [53] 

k-€ Model 

X R 

x R = 5.8 (Low Re flow) 
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Ilegbusi and Spalding [51] with their modified k-co model obtain a value of = 7.2 H. Their 

conclusion, however, is that both k-e and k-co model (modified) predict the flow over a backward facing 
step equally well. The difference between the predicted and measured peak shear stress across the channel 
ranged from -17.5 percent to -1 percent for the k-e model as compared to -4 percent to +21 percent 
for the k-co model. Static pressure is well predicted by both the models. 

Launder [56] discusses the effects of mesh size, false diffusion and inlet boundary conditions on 
the flow prediction. These are considered in Section VI. 

There arises a question about the prediction of the flow over a backward facing step. Why is it 
that the flow in a pipe expansion is predicted well (Section 5.2.1) but not the asymmetrical channel 
expansion? One is reminded of the problem of the prediction of plane and axisymmetric jet flows 
[57,58], where the spreading rate of a plane jet is correctly predicted while that of an axisymmetric jet 
is overpredicted using the standard k-e model. It is not exactly clear what causes this difference in 
predictions between the two confined flow geometries. An attempt is made to explain the behavior of 
the k-e model in the two cases in Section 5.2.1. 


5.1.2 Symmetric Expansion (Double Step) Flow 

Gosman et al. [59] studied the plane symmetric sudden expansion (sometimes called a double 
step) flow using k-e model. The regions of recirculation on both sides of the symmetric plane were 
identical for low expansion ratios. For large expansion ratios (>1.5), they observed that the two recir- 
culation regions interact to produce asymmetric regions of recirculation. At some conditions, they 
observed two regions of recirculation on one side and one on the other side of the symmetric plane. 
The numerical results for the range of expansion ratios where symmetry was observed showed good 
agreement with experiments [60]. Though one would expect that the stabilizing effect of the top wall 
(instead of a symmetry line) to produce a shorter recirculation region in the case of asymmetric (single 
step) expansion flow as compared to double step flow, it has not been precisely measured in the 
experiments. 


5.1.3 Flow Over a Square Obstacle 

Durst and Rastogi [61] investigated theoretically as well as experimentally (using LDA) the flow 
over a square obstacle in a two-dimensional channel. They used a k-e model and a three equation k, e 
and uv model. They obtained a recirculation region behind the obstacle of size x^ = 7.5 H. The calcu- 
lations showed reattachment on the obstacle itself (Fig. 8) which was not observed in the experiments. 
This they attribute to three-dimensional effects. Velocity profiles upstream were faithfully reproduced, 
while over and behind the obstruction the agreement was poor. The redevelopment after the reattach- 
ment was also slow. Use of the three equation model improves the calculated mean velocity profiles, 
but not significantly. The kinetic energy profiles show only qualitative agreement. The calculated length 
scales in the separated region and downstream of separation point indicate large scale eddies. They 
observed that the k-e model has to be modified to obtain more accurate calculations in the separated 
region. 
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Figure 8. Flow over a square obstacle: predicted and measured flow pattern [61]. 

5.2 Axisymmetric Flows 

5.2.1 Flow in a Sudden Pipe Expansion 

Flow in a pipe expansion is one of the basic problems encountered in common engineering flow 
systems. The separated, reattached and recirculating flow downstream of the expansion greatly influences 
the transport of momentum and heat. Increased levels of heat and mass transfer occur. This increase is 
attributed to increase in levels of turbulent kinetic energy in the stream. The high kinetic energy comes 
about due to the large shear stress in the separated layer and the associated low rate of turbulent energy 
dissipation. Because of the intrinsic practical interest several experimental and theoretical studies have 
been carried out on this flow configuration. This geometry also forms the basic configuration for combus- 
tor studies (dump combustors), with swirl introduced before the expansion. In earlier computations of 
pipe expansion flow, a uniform length scale over the entire separated region was assumed [11]. Only 
with the introduction of a differential equation for calculating the variation of the length scale in the 
separated region, computations can predict the desired flow properties (Section 5.1). 

As in the backward facing step flow, an important parameter for comparison of pipe expansion 
flow predictions, is the reattachment length x^. The location of the reattachment point also determines 

the location of maximum heat/mass transfer coefficient. Figure 9 shows the predicted separation lengths 
by different computations. Table 2 lists the turbulence model used by different authors and some 
observations. The experimentally observed value of x^ is 8.5 to 9 H for this geometry [64,68,69]. From 

the figure one notices the following: In contrast to the predictions of backward facing step flow, the 

standard k-e model predicts the reattachment length within the experimental uncertainty. This is 
interesting since it has some implications on the set of model constants used. One can discern the follow- 
ing prediction trends when one uses the same model constants appropriate for the boundary layer flows. 

The flow in a sudden pipe expansion is predicted fairly well while the flow over a backward 
facing step is underpredicted by about 20 percent, using the standard k-e model. It is not clear whether 
this is related to the problem of the prediction of plane and axisymmetric unconfined jet flows [57,58], 
However, the underprediction of the flow over a backward facing step flow may be related to the under- 
prediction of the recirculation region in the bluff-body stabilized flows [70], Habib and Whitelaw 
[71,72] observe that the k-e model connects the dissipation rate/production strongly to the mean flow 
field. In regions where there is more than one major component of velocity gradient tensor extra strain 
rates are added to the generation terms. Then, the Reynolds stresses change by values that are much 
larger than the direct effect of these strains and this is not represented by the effective viscosity model 
[32]. The response of the k-e model to changes in mean flow field is too fast. The stabilizing effect of 
the top wall in the flow over a backward facing step may not be correctly represented in the k-e model. 
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Figure 9. Axisymmetric expansion; predicted separation length, x^. 

Habib and Whitelaw [71] find that in the confined coaxial jet predictions, with the increase in 
the velocity ratio of annulus to pipe greater than unity, the underprediction of the separation region 
increases. The top wall in the backward step flow may have similar effects. 

Reverting to Table 2, the modified k-t models of References 66, 67, 68, and 28 do not appre- 
ciably improve the reattachment length. Chieng and Launaer [28] state that the near wall model does 
not change the reattachment region but only improves the heat transfer prediction. Tire reduced value 
of xr in these cases may be partly due to the low Reynolds number (R e [)) of these computations. 

Reference 64 gives the variation of x^ with R e p). 

To conclude this section, it is observed that no standard modification of the k-e model which 
will respond correctly to the axisymmetric as well as asymmetric plane expansion flow has been arrived 
at. It appears that the algebraic stress model would be a better alternative. 
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TABLE 2. FLOW IN A SUDDEN PIPE EXPANSION 


i 
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5.2.2 Diffuser Flows 


The axisymmetric diffuser flows can be divided into two groups: (i) small angle diffusers (no 
separation) with expansion of total angles equal to 5 to 10 deg, represent the aerodynamic diffusers; 
(ii) large angle diffusers (with separation) as found in “dump” combustors (90 deg half angle represents 
the sudden expansion). The computation of group (i) flows is generally through the solution of parabolic 
type equations and will not be considered here. In the second class only two recent contributions will be 
discussed. 

Habib and Whitelaw [73] computed the axisymmetric recirculating flows in wide angle diffusers 
and compared this with their experimental results. They considered 20 and 40 deg half-angle geometries. 
Tire diameter expansion ratio was 2.0. They computed the flow with an orthogonal-curvilinear (body 
fitted) coordinate system. The standard k-e model was employed. Mean axial velocity and radial velocity 
profiles and kinetic energy contours were presented and compared with experiments. 

They find that the location of the maximum kinetic energy is correctly predicted but its value 
is underpredicted by about 30 percent. They attribute this to the incorrect representation of the source 
terms in the transport equation for e and partly to the extra strain terms in the calculations of Reynolds 
stress terms and the rate of dissipation. The mean axial velocity and radial velocity are well predicted. 
Rhode et al. [74] report the computed results for expansion angles of 45, 70 and 90 deg, for the 
diameter expansion of 2. They use a non-uniform grid system and the standard k-e model. The predicted 
mean velocity profiles show good agreement with the experiments of [69], The computed flow pattern 
agrees with their flow visualization pictures. They characterize the size and shape of the recirculation 
region as a function of the angle of the sloping wall. They, however, did not make any comparison of 
the turbulence quantities such as the Reynolds stress. 


5.3 Curved Flows 

This section discusses the predictions of flows with streamline curvature caused by the surface 
curvature and swirl. 


5.3.1 Flows with Surface Curvature 

Pourahmadi and Humphrey [34], using their modifications for the value of C^ discussed in 

Section 4.3, obtained very good predictions for the fully developed flow in a highly curved two- 
dimensional duct. 

Figure 10 shows the streamwise variation of friction factor at the inner and outer walls of 
strongly curved channel flow. The use of standard C^ in k-e overpredicts the innerwall Cf while severely 

underpredicting the outer wall Cf. The modified model predicts the inner wall Cf very well and greatly 
improves the Cf at the outer wall. These predictions are in agreement with the experiments ofHonaini 

et al. [75], The prediction of the mean velocity profiles is satisfactory even with standard k-e models. 
Launder et al. [35] predicted the curved duct flow corresponding to the experiments of Lillis and 
Joubert [76], with the modification of lengthscale incorporating the gradient Richardson number. They 
obtain good agreement for the mean velocity and the shear stress at the inner wall. 

Towne [77] has computed the turbulent flow through circular and square cross sectioned S-ducts 
using three-dimensional parabolic equations with marching procedure. He obtains good agreement of the 
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mean velocity profile with experiment, but no details of the shear stress variation is reported. He used 
an algebraic mixing length turbulence model. Typical streamwise velocity profiles in the pipe are shown 
in Figure 11. It is seen that a fine mesh (50 x 50 x 80) is needed to produce a meaningful prediction of 
the mean velocity profile variation. 



Figure 10. Streamwise variation of friction factor at the inner and outer walls 
of strongly curved channel flow [341. 



EXPERIMENT 
ANALYSIS, 50 x 50 x 80 
ANALYSIS, 25 x 25 x 80 


Figure 11. Computed and measured streamwise velocity' profiles in symmetry plane 
for turbulent flow in 22.5 to 22.5 circular S-duct [77], 
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Mukeijee et al. [78] use the k-e model to compute the flow in a highly curved 180 deg turnaround 
axisymmetric duct configuration of the SSME. The curvature effect has not been explicitly included in 
their turbulence model. Because of this, it is difficult to estimate the accuracy of the solution, though it 
appears to predict the trend qualitatively [79]. No detailed predictions of the turbulent flow in a highly 
curved channel has been reported in the literature. Agarwal et al. [80] using LDA obtained detailed 
laminar flow characteristics in a 180 deg curved tube, of R c /r = 7 and 20. This flow was successfully 
simulated using a three-dimensional code with fine mesh (50 x 50 x 226) by Towne [77]. No detailed 
experimental results are available for the 180 deg bend of R c /r se 1 as found in SSME turnaround ducts. 


5.3.2 Flows with Swirl 

Swirling flows result from the application of a tangential velocity component imparted to the flow 
by swirl vanes, tangential entry, etc. Swirl flow occurs in a variety of applications, the most 
extensively studied one being the combustion aerodynamics. Nuttal [81] was the first one to observe 
the reverse flow characteristics of a swirling flow in a circular pipe. For Reynolds numbers in the range 

1 x 10 4 to 3 x 10 4 , three types of flow patterns are observed as shown in Figure 12. Transition from 
one-celled to two and finally a three-celled vortex structure occurs, as the swirl component of the 
velocity increases. This would be a difficult problem for modeling. A compendium of experimental and 
computational studies of swirl flows can be found in the recent book, “Swirl Flows,” by Gupta et al. 
[82]. The review by Jones and Whitelaw [83] provides a good account of the prediction methods for 
reacting flows. The following sections discuss predictions of two nonreacting swirl flow configurations, 
namely the swirl flow in a pipe expansion and the swirled confined coaxial jets. 



ONE CELL 
LOW SWIRL 




TWO CELLS THREE CELLS 

MEDIUM SWIRL HIGH SWIRL 


Figure 12. Flow reversal observed by Nuttal for swirling flow in a circular pipe [81], 


5.3.2. 1 Swirl Flow in a Pipe Expansion 

This is the basic form of the dump combustor geometry. With the introduction ol swirl, 
a central recirculation zone (CTRZ) is formed, in addition to the corner recirculation zone (CRZ) 
[74]. The recirculation zones are important in the design because most combustion occurs in 
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and near these regions. The nonswirling pipe expansion has only the CRZ and with the swirl imparted at 
the inlet, a central recirculation region appears and the size of CRZ decreases. The comer recirculation 
zone is reduced in the axial extent until it disappears for a swirl vane angle of 70 deg (Fig. 13). These 
predictions were obtained using the TEACH code with k-e turbulence model. Novick et al. [84,85] obtained 
useful flow predictions for actual combustor geometries. One example is shown in Figure 14. With no swirl 
the CRZ fills the outer cavity. As swirl is introduced a CTRZ appears in conjunction with a decrease in the 
size of CRZ (see Fig. 14 for 45 deg swirl angle). Further increase in swirl angle results in continued enlarge- 
ment of the CTRZ and a decrease in the size of the CRZ. At 75 deg the CRZ disappears. No direct com- 
parison with experiments has been made. 




^ ; 1 : 

0.00 1.00 2.00 3.00 

(b) <t> = 45° 



AXIAL POSITION x/D 

Figure 13. Predicted streamline plots with wall expansion angle a = 90 deg 
and various swirl vane angles 0. [74]. 
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5. 3. 2. 2 Confined Coaxial Jets 


This basic configuration is being studied extensively both experimentally and theoretically because 
of its importance in combustor designs. Efforts have been made to predict the flow patterns for a range 
of swirl numbers. At low swirl numbers S < 0.5, there is only the comer recirculation zone. For swirl 
numbers S > 0.5 an additional central recirculation zone is observed. The swirl number is defined as 


R 

f p Ur 2 W dr 



R / p (U 2 - 0.5 W 2 ) r dr 
0 


(S Q is the swirl number of the outer jet, S- is the swirl number of the inner jet, and a is the velocity ratio of 
annulus to pipe). 

Table 3 summarizes the recent numerical studies on swirling jets. As it is obvious from the table, 
all the studies use the k-e turbulence model. Detailed analysis of the predictions of the k-e model for 
swirling jet is presented by Habib and Whitelaw [72]. The recirculation length is underpredicted by about 
20 percent for a = 1. With increase in a, the prediction gets worse. This is attributed to the incorrect 
representation of the diffusion process in the k-e model. For the flow with swirl, the model is incapable 
of predicting the velocity minimum obtained at the axis in the experiments. Figure 1 5 shows the velocity 
profiles for swirl and no swirl cases. This velocity minimum on the central line has not been predicted 
correctly in any of the studies reported. They also observe that the intermediate swirl number of 0.23 
provides the flow which is more difficult to present by calculation methods since the near recirculation 
region is located away from the solid surfaces. With no swirl there is only the CRZ. With the swirl 
number of 0.5, the CTRZ is tied to the exit geometry and the downstream flow is easily represented. 
Srinivasan and Mongia [41] modified the k-e model to include the Richardson number dependence 
through the constant LA of the model. They use the gradient Richardson number, and divide this into 

two: swirl Richardson number and curvature Richardson number quite arbitrarily. However, they find 
that at high swirl numbers, the curvature Richardson number has no influence on the flow’ but swirl 
Richardson number is the controlling one. They could not obtain convergent results for coswirl case. 
They modified their program to include radial pressure variations, which enabled them to obtain coswirl solu- 
tions. Coswirl solutions showed a central recirculation zone not observed in experiments. Ramos [86] 
claims that with the use of suitable initial condition, the standard k-e model predicts the coswirl flow. 
His claim is not substantiated by providing the exact initial condition used and the justification for the 
conditions used. Hence, it is not of much use in understanding the prediction method. Novick et al 
[84,85] study several variations of the combustor configurations. They characterize the size and shape 
of the recirculation zones with swirl number, expansion ratio and the geometry. These results help one 
to get some idea of the predictive capability of the k-e turbulence model for these configurations. 

In general, in the prediction of confined vortex flows, the k-e model performs rather poorly [88]. 
Figure 16 shows the predicted and measured tangential velocity profile in a vortex tube. It is seen that 
the standard k-e model does not even predict the qualitative features of the flow while the ASM pre- 
dictions agree well with experimental results [89]. Dixon et al. [90] employed the k-e model for pre- 
dicting the swirled coaxial jets and found that the model underpredicted the extent of the central recircu- 
lation region by about 20 percent. They also found that the prediction of the maximum recirculated mass 
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TABLE 3. CONFINED COAXIAL JETS 
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flow rate in the central recirculation zone was still worse; it was half that of the measured one. Thus, 
the use of the standard k-e model to highly swirling flows may not produce reliable results [91 ] and the use 
of ASM should be recommended. 

The accuracy of the numerical schemes in predicting the swirl flows and the effect of inlet 
boundary conditions are briefly discussed in Section VI. 




Figure 15. Measured and calculated distributions of the mean velocity, o Measurements, 

— calculations, (a) S = 0.0. (b) S = 0.23 [72] . 



Figure 16. Predicted and experimental tangential velocity profiles in the vortex tube. 
- ASM, — k-e model, a experimental. [89] 


35 




5.4 Three-Dimensional Computations 
5.4.1 Non circular Duct Flows 

Turbulent flow in noncircular ducts is characterized by the everpresent secondary motions. The 
magnitude of the secondary motion depends on the geometry and is strongly influenced by the curvature. 
Presence of secondary flows makes the problem more difficult. One has to solve for the three 
components of velocity and additional equations for Reynolds stresses. To make the problem tractable, 
most studies on non circular passages concentrate on fully-developed flow in uniform cross-section. 
(Developing flow was a test problem in the 1981-Stanford conference and will be presented in this 
section.) 

The cause of secondary motion has been investigated and explained. Recently, Speziale [921 
has proven mathematically that the secondary flows in noncircular ducts result from a nonzero difference 
in the normal Reynolds stresses on planes perpendicular to the axial flow direction. He also shows that 
the widely used k-e turbulence model has no built-in mechanism for the development of secondary flow, 
while the second order closure models do. 

Launder and Ying [93] were the first to carry out a detailed numerical prediction of the turbu- 
lent flow in a square duct. They employed the simplified form of the Reynolds stress model of Hanjalic 
and Launder [94] which advocates the solution of seven coupled nonlinear differential equations in 
addition to those of mean flow quantities (Section 3.3). Launder and Ying found that for fully developed 
flows the convection and diffusion terms in the stress transport equations may be neglected. This 
important simplification results in the removal of all differential coefficients of ujuj and hence algebraic 

in the latter (the reasoning putforth in ASM, Section 3.5). They obtain the following form of equations: 



They solve the kinetic energy equation in addition to the three momentum equations and use the 
algebraic relations for the stresses (49) and the length scale. 

Later investigations employed a differential equation for e. This method of solving five differential 
and algebraic equations for the stresses has been successfully employed for solving the flow in rectangular, 
triangular and elliptic ducts, and tube assemblies by Gosman and Rapley [95], They use a boundary fitted 
(orthogonal curvilinear) coordinate system for the noncircular ducts. Comparison of the ASM equations 
with laminar contributions reveals that axial stresses UjU^ and uyu^ depend on turbulent viscosity and 

axial strain rate in a Newtonian fashion; but the cross-planar stress equations are distinctly different from 
the Newtonian relations depending on the axial as opposed to cross-planar strain rates. It is this feature, 
as explained by Speziale, that is responsible for the turbulence-driven secondary motion. 

Perhaps, the most exhaustive study of fully-developed turbulent flow in noncircular ducts is that 
of Gosman and Rapley [95], They studied the flow in square, rectangular, triangular, and elliptic ducts 
and tube assembly passages and compared the accuracies of the predictions by different turbulence 
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i 

I 

| models. Detailed comparison of the mean velocity, secondary velocities, friction factor and Reynolds 

I stress components with those of experiments [96] and other numerical predictions were presented for 

the square duct. They compared the prediction of the algebraic stress model, Reynolds stress model [97], 
and modified Reynolds stress model [98]. The modification in Noat et al. [98] is only the finite differ- 
ence method of solution. By employing a local block inversion method, they solved the flow and stress 

equations simultaneously rather than sequentially as in other methods. Figures 17 and 18 show the com- 
parison of the prediction by the three models with the experiment [99], The three models predict the 
axial (Fig. 17) and secondary velocities (Fig. 18) equally well. However, the Reynolds stress components 
(Fig. 19) are not predicted well. This may be partly due to the expected underprediction of the normal 



Figure 18. Secondary velocity profiles (Uj/U^*) in a square duct, Re - 2.15 x 10^ [95] . 
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o U 3 \ U-,' 2 , x Uf U 3 EXPERIMENT 
BRUNDRETT AND BAINES [99] 

— PREDICTION, GOSMAN AND RAPLEY [95] 


Figure 19. Reynolds stress profiles in a square duct, Re = 8.3 x 10^ [95]. 


5.4.2 Developing Flow in a Square Duct 

This was a test problem at the 1981 -Stanford conference on complex turbulent flows. Only three 
groups attempted this problem. Table 4 shows the models used by them and comments on the predicted 
results. All three groups use the ASM-based method. References 46 and 100 use a differential equation 
for e, while Reference 101 uses an algebraic equation. As indicated in the table, the model predictions 
are very poor. This means that the algebraic stress model is not good for developing flow where local 
equilibrium assumption is not valid. Perhaps only the Reynolds stress model is capable of predicting the 
developing flow in a square duct. 
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TABLE 4. DEVELOPING FLOW IN A SQUARE DUCT 
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VI. NUMERICAL METHODS 


The prediction methods use various numerical schemes and solution procedures whose perform- 
ance in different regions of the same flow and in different flow situations are not well established. 
Two related problems which influence the predicted solutions are briefly discussed below. 


6.1 False Diffusion 

Most methods employ an upwind difference scheme which is known to introduce false diffusion 
errors. The false diffusion coefficient at any computational cell can be computed from the formula 
[102,103] 

p U AxAy sin 29 

r f = — ^ 5 — (50) 

4 (Ay sur 9 + Ax cos'’ 6) 


where Ax and Ay are the cell dimensions and 9 is the angle made by the resultant velocity at the grid 
lines. Figure 20 shows, for a hybrid central-upwind scheme with a fine grid (42 x 42) which produces a mesh 
independent solution, the variation of maximum normalized false diffusion coefficient along the duct of 
backstep flow (Section 5.1). The false diffusion coefficient shown is the value normalized with local value of 
turbulent viscosity and multiplied by the ratio of the local-to-maximum shear stress at any station [56] . The 
maximum value of the factor is about 0.2 and the estimated error caused by this false diffusion in the pre- 
diction of the reattachment length is about 0.3H. For example Syed et al. [53] , by using an accurate scheme 
(bounded Skewed Upwind Scheme), obtain an increase of 0 . 35 H in x^ over hybrid scheme, for a fine 

mesh. In the presence of numerical diffusion in the recirculation region, no meaningful comparison of 
two turbulence models can be made. On a moderate grid, upwind differencing scheme may be 
insensitive to changes in turbulent viscosity. If one is unable to accommodate a fine enough mesh, then 
the alternative is to use a higher order accurate scheme. McGuirk et al. [104] use a selective mesh refine- 
ment procedure which reduces the number of mesh points required but maintaining a higher accuracy. 
They identify the regions of “false diffusion” present in the solution, with a given grid, by estimating 
the magnitude of the local truncation error in the converged solution of the discretized equations. This 
identifies the local regions where the Peclet number is high. They change the difference scheme depend- 
ing on the value of the Peclet number: if |Pe| < 2, they use central differencing scheme and if 
|Pe| > 2, they use upwind difference scheme. This guarantees numerical stability and produces accurate 
results. This procedure may be useful in comparing two different turbulence models. But, it is too involved 
and expensive to use as a general solution procedure. 

Numerical differencing schemes that reduce the false diffusion have been developed [105,106]. 
Leschziner and Rodi [37] evaluated the performance of three difference schemes - hybrid central/ 
upwind differencing scheme (CUDS), hybrid central/skewed-upwind differencing scheme (CSUDS), and 
quadratic upwind weighted differencing scheme (QUDS) — on annular and twin parallel jet flows and 
found that the later two (CSUDS and QUDS) reduce the false diffusion errors significantly. Syed and 
Chiappetta [53] have tested several schemes for accuracy in predicting the swirl flows. They found that 
the bounded skewed upwind scheme (BSUDS) produced more accurate results than the hybrid (CUDS) 
scheme. 
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Figure 20. Maximum relative level of false diffusion coefficient at any section [56]. 


6.2 Inlet Boundary Conditions 


Inlet conditions are usually considered of minor importance while evaluating the performance of 
turbulence models against experimental data. Assumed boundary conditions can lead to wrong conclu- 
sions about the performance of the model. In particular, the inlet profiles of k and e can have a 
significant effect on the flow downstream. The profiles of k may be estimated from the experimental 
results. But no satisfactory method of evaluating inlet e profile is available. In the flow over a backward 
step, for example, if one starts with a larger initial length scale (larger viscosities), the prediction is 
expected to give a shorter reattachment length. But the predictions do not show such a trend. 


The inlet conditions become more “controlling” when the inlet flow has a swirl component. 
Sturgess et al. [87] studied in detail the effect of inlet boundary conditions in swirl flows. They found 
that the computations were sensitive to the inlet conditions, and might affect the solution in the entire flow 
domain. Many times, by manipulating the inlet profiles, a better fit to the experimental data can be 
obtained. For example, Ramos [86] claims that k-e model predicts a recirculation zone for both co- and 
counter-swirl conditions, if suitable initial conditions are used. Leschziner and Rodi [107] find that the 
inlet conditions for k and e play as crucial a role in achieving predictive accuracy as turbulence modeling 
details. Thus, in the absence of experimental inlet profile data, careful and “judicious” choice of the 
initial condition estimation is essential for a physically meaningful prediction. 
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Typical Inlet Conditions Employed in Predictions 



k in 

e in 

(a) Flows Without Swirl 

Qin [62] 

0.003 x U in 2 

0.09 x k in 3 / 2 /(0.03 x D/2) 

(b) Flows With Swirl 

Habib and Whitelaw [72] 

k = 1/2 [u' 2 + v' 2 + w' 2 ] 

0.09 x k in 3 / 2 (0.03 x R) 


u' - experiment 

(R is the radius of pipe or 


v' = w' = 0.6 u 

annulus) 


Sturgess et al [87] studied the effect of inlet boundary conditions on the prediction of swirl 
flows. They found that the numerical solutions were very sensitive to the inlet boundary conditions and 
the overall accuracy of the prediction depended on the grid system. Syed and Sturgess [108] describe the 
present status of the predictive capability for the swirling recirculating nonreactive flows (Fig. 21). To 
this should be added the inability to predict the coswirl flow, the central line minimum, and the radial 
variations of the mean velocity as discussed in Section 5.2.3. Detailed experimental data for coaxial jets 
with and without swirl can be found in References 109 and 110, for evaluating turbulence model pre- 
dictions. 


REGION I 


REGION II 


REGION III 


INLET REGION: 

ACCURACY OF 
PREDICTIONS 
DETERMINED BY 
ACCURACY OF 
INLET PROFILE 


RECIRCULATION 
REGION: 
STRENGTH = 
UNDERPREDICTED 
BY 20 - 25% 

SIZE = 

UNDERPREDICTED 
BY 15 - 20% 


PARABOLIC 

REGION: 

ACCURACY BETTER 
THAN 85 - 90% 


PROFILE = 

QUALITATIVELY 

PREDICTED 



Figure 21. Summary of the accuracy of predictions in the various regions of flow with swirl [108]. 
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VII. COUNTERGRADIENT TRANSPORT AND TURBULENCE MODELING 


The eddy viscosity type models assume that the gradient transport hypothesis holds in all the 
regions of the flow. However, it is known that small regions of countergradient transport or negative 
turbulent energy production exist in many flow situations such as asymmetric channel flow, wall jets, 
wake flows, perturbed shear layers, etc. Beguier et al. [Ill] discuss a variety of flow situations which 
contain regions where the product -uv 9U/3y becomes negative. Townsend [112] was the first to 
recognize and explain the small region of countergradient transport of momentum in a plane wake. 
In this region, transport of momentum consists of two parts: (a) gradient type diffusion by small 

scales and (b) convection by large scale motions of eddies comparable in size to the characteristic 
length scale of the problem. Hinze [113] extended this concept to asymmetric flow. To account for the 
contribution due to the bulk transport by the large eddy, one adds a second term in the expression for 
Reynolds stress (see, for example Ref. 3, page 372). The additional term is proportional to either 
3( v 2)l/2/gy or 0(q2)l/2/gy w h ere (q2)/2 is the kinetic energy. 

Beguir et al. [Ill] analyzed different flow situations in which such regions of countergradient 

transport occurs. They find that a direct relationship between the sign of the gradient of (q^/2 
and the displacement, e, of the minimum exists. The displacement e is defined as 


e y (uv=0)~ y (9U/9y=0) 


— 1 0 3 (q 2 )'P 

v* u 3y 


(51) 


where C is a constant (= 0.1), 1 Q is a characteristic length scale, such as the width of the channel. The 

product e • 9(q^)^/9y was negative for the majority of the flows considered, where the flow asymme- 
try or the perturbation was small. However, when the asymmetry, such as the interaction of wake behind 
cylinders of unequal diameters, or the perturbation was strong, the product was positive. Their flow 
visualization of the interacting wake behind the cylinders showed that in that case, merger (or pairing) 
of the vortices of the same sense of rotation occurred, similar to that found in shear layers [114,115]. 
This suggests that the small asymmetry or perturbation results only in a change in the orientation of the 
generally elliptic vortex structure (against the gradient 9U/9y instead of along the gradient) and results 
in only minor changes in turbulence quantities. This results in a negative sign of the product 

e • 9(q^)^^/9y. The strongly asymmetric or perturbed flow condition results in merger of large scale 
eddies, in which case the variation in turbulence quantities is too complex. This situation cannot possibly 
be described by simply adding a term to the Reynolds stress equation. An example of the changes in tur- 
bulence quantities in a perturbed axisymmetric mixing layer [116] is shown in Figure 22. It can be seen 
that for the strong perturbation, in which localized vortex merger occurs, the shapes of the Reynolds 
stress and the longitudinal velocity fluctuation profiles are different from those of unperturbed or of 
small perturbation and are too complex. Oster and Wygnanski [117] present a detailed study of the 
variations of the turbulence quantities in a perturbed plane mixing layer. Employing Reynolds stress 
model one can predict the small asymmetric or perturbation case satisfactorily [23] or even by a 
simplified form of the Reynolds stress model [94], However, it is not known if the Reynolds stress 
model can predict the strongly asymmetric or perturbed case. 
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Figure 22. Perturbed mixing layer characteristics [116]. 

For the transport of any scalar quantity, 0 in a turbulent flow, one can define an equation similar 
to that of (51). 


= „ _ C 0 , 3 9(q 2 ) 1//2 

H y0v=o y(90/3y=o) 0 


Beguier et al. [Ill] considered this for the case of asymmetric mean temperature ( 6 ) profile. 

Again for small perturbations the product e^ • d(q 2 ) , / 2 /3y was negative and for a strong interaction 
case it was positive. 






The countergradient diffusion in flows involving combustion has been studied by Libby and his 
co-workers [118,119]. They avoid the gradient transport assumption and develop a second order closure 
model. They attribute this diffusion to the differential effect of the mean pressure gradient on cold 
reactants and hot products. Jones and Whitelaw [83] also observe that the countergradient diffusion is 
due to the preferential influence of the mean press ure g radient on low- and high-density gases which 
manifest through the terms like p'uj" 0p/dx- and p'$" 3p/9xj in the exact Favre-averaged (density 

weighted) [120] Reynolds stress and turbulent scalar flux transport equations. 

Spalding [121] calls the countergradient diffusion as “pressure gradient” diffusion, because of its 
connection to the pressure gradient. He observes that the turbulence modelers have been concerned exclu- 
sively with Kelvin-Helmholtz (shear generated vortex roll-up and pairing) instability and little attention 
has been given to Rayleigh-Taylor instability. His two-fluid model takes care of the later and is expected 
to describe the so-called countergradient or pressure gradient diffusion process. However, it is not clear 
at this point, if it can describe the shear generated variations of turbulence quantities in small and 
strongly perturbed/asymmetric flow discussed above, as well. 


VIII. COMPUTATIONS OF RECIRCULATING FLOWS WITH PHOENICS CODE 


PHOENICS is a general purpose flow simulation code developed by Prof. Spalding and his asso- 
ciates at the Imperial College, London. The acronym PHOENICS stands for Parabolic Hyperbolic or 
Elliptic Numerical Integration Code Series. The code is now available commercially through the CHAM 
of North America in the U.S. and the CHAM of the United Kingdom in the U.K. The general features 
of the program are discussed by Spalding in Reference 122. The basic concepts and solution procedures 
employed in the code can be found in Reference 123. Here we present the computational results 
obtained with the k-e model using the PHOENICS code for two flow configurations: (i) the flow over 
a backward facing step and (ii) the flow over a square obstacle in a two-dimensional channel. Only a 
brief discussion on the results is included. The sensitivity of the solution to the inlet boundary condi- 
tion and the difference in the k-e model predictions for the axisymmetric and asymmetric channel 
expansions will be discussed in a forthcoming report. 


8.1 Flow Over a Backward Facing Step 

This is the benchmark problem discussed in Section 5.1.1. The layout of the grid for the present 
computation is shown in Figure 23. The grid used is an adjusted (more grid lines near the walls) 42 x 42 
grid to produce a grid independent solution. The expansion ratio is 1-5. The flow pattern with the recir- 
culating region behind the step is shown in Figure 24. The length of the recirculation region x^ is 

equal to 5.8H, about the same as obtained in other k-e model predictions. Depending on the ratio of 
the mesh sizes in the x and y directions (Ax/Ay), sometimes the solution exhibits a long tail for the 
recirculation region [62], Figure 25 shows the contours of the kinetic energy in the flow domain. A large 
variation of the kinetic energy is observed in the separated shear layer and in the wall layer downstream of 
the step (the contours in the core region have very low levels of k as would be expected). The contours of 
e, the dissipation rate of k, are shown in Figure 26. Here again contours in the core region have very low 
levels of e. The maximum variation of e occurs near the convex corner of the step. Figure 27 shows the 
static pressure contours in the flow field. As would be expected, the maximum static pressure variation 
occurs across the recirculating region. Thus, the predicted flow parameters show qualitative agreement with 
the experimental results. But, quantitatively, the prediction shows poor agreement with the measured ones 
both in the recirculating and redeveloping regions. 
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Figure 26. e-contours. 
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Figure 27. Static pressure contours. 


8.2 Flow Over a Square Obstacle in a Two-Dimensional Channel 

The flow geometry and the flow conditions are the same as in Reference 61 and discussed in 
Section 5.1.3. The grid layout for this problem is shown in Figure 28. The grid lines are uniformly 
spaced in the y-direction and adjusted in the x-direction. The grid employed is 40 x 47, which produces 
a nearly grid independent solution. There are four grid lines in the width of the obstacle. The flow 
pattern obtained is shown in Figure 29. The flow separates at the upstream comer of the obstacle and 
there appears to be a small recirculating region (indicated by the sharp turn of this streamline and the 
static pressure contours, Fig. 32) on the obstacle. Durst and Rastogi [61] found this separating streamline 
to reattach on the obstacle itself (Fig. 8). However, their experimental results show no reattachment on 
the obstacle. The length of the separated region is x^ = 6.6 H in the present prediction as compared to 

~ 7-5 H observed in the measurement. Considering the nature of the flow, this flow is slightly better 
predicted than the flow over a backward facing step. The reattachment length is underpredicted by only 
about 12 percent, while for the backstep flow is underpredicted by about 20 percent. The kinetic 

energy contours are shown in Figure 30. Because of the accelerating and decelerating nature of the flow 
no core region (downstream of the obstacle) in which the contour levels are very low, is observed, as in 
the backstep flow. Figure 31 shows the contours of the dissipation rate e in the computational domain. 
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The static pressure contours shown in Figure 32 indicate large variations of pressure in the recirculation 
region behind the obstruction and the region over the obstruction. Thus the qualitative features of the flow 
are in good agreement with experiments. Further aspects of the predictions of these flows will be reported 
in the forthcoming report. 



Figure 28. Calculation domain and grid layout for flow over a square obstacle. 



Figure 29. Streamline pattern. 
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Figure 30. 


Kinetic energy contours. 
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Figure 31. e-contours. 
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IX. OTHER APPROACHES TO TURBULENT FLOW COMPUTATION 


In this section, brief descriptions of three other approaches to turbulent flow computation, 
namely the vortex method, the large eddy simulation and the direct simulation are given. These methods 
are sometimes called higher level simulations, since they incorporate more flow physics than is possible 
with turbulence models. They are still more of a research tool than a scheme for engineering calcula- 
tions. The computational memory and time requirements of these methods are typically high. Hence, the 
time for their implementation to aid the engineering flow analysis is yet to come. 


9.1 Vortex Method 

With the recognition of the existence of large-scale organized structures at least in the initial 
regions of mixing layers, jets and wakes, vortex method has received renewed emphasis in recent years. 
It has been used extensively to study the dynamics of large scale structures, since these large scale 
features arise from the properties of a rotational but inviscid flow. The idea is to represent turbulence 
as a superposition of interacting vortices. Saffman and Baker [124] provide an excellent account of 
vortex interaction studies. Leonard [125] reviews the concepts and methods of vortex flow simulation, 
including those for three-dimensional flows. Aref in his recent review [126] on the solution of two- 
dimensional Euler equations, discusses the use of vortex methods in the studies of transition to turbu- 
lence and the onset of chaos. 

The widely used vortex method is the point vortex method. In this approach, the vorticity 
originally confined in a thin layer is concentrated further into a finite number of point vortices. In other 
words, the piecewise continuous distribution of vorticity co is replaced by a sum of N 5 functions. 


co(X,t)=]T r n 5[X-X n (t)] (53) 

n=T 

where 5 is the Dirac delta function, X n = (x n ,y n ) is the location of the nth vortex of circulation T n . 
The motion of vortices is followed by integrating the system of ordinary differential equations 


dX„ 


dt 


= U(X n ,t) 


(54) 
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The solution of equations (53) and (54) represents the solution of the two-dimensional Euler equations: 


dco 

— + (U • V ) co = 0 
dt 

j 

(55) 

V ^ ^ = -co ) 


(56) 

d\p 

u = — . U = 

dp 

(57) 

t 

X 

QJ> 

<< 

dx 


No precise answer can be given for the question of the number N of the vortices required for a simula- 
tion. For a given N and t, confidence in the method is based on the requirement of the simulation and 
is often limited by the computer memory and time requirements. 

Experimental studies on perturbed mixing layers and jets suggest that the sensitivity of the 
normally turbulent flows to external forcing is related to the rotational in viscid behavior of these flows. 
Hence, vortex methods have been used to study these perturbed flows. Acton [127] studied the charac- 
teristics of the “preferred mode” in circular jets. Vortex method has also been used to study the turbu- 
lence suppression in perturbed shear layers, and vortex roll-up, pairing, etc., in perturbed plane jet flows 
[128,129]. These studies depict clearly the dynamic features of the large scale vortical structures 
observed in the experiments. The vorticity and Reynolds stress contours during pairing of vortex struc- 
tures in a perturbed mixing layer obtained using the vortex-in-cell method are shown in Figure 33. This 
corresponds to the orientation of pairing structures, for which a region of countergradient transport 
will result in the time mean field. 

In Chorin’s method for viscous flows [130], the inviscid outer flow is solved via discrete vortices 
and viscous effects are incorporated via random walks for each vortex. Recently Chorin and Iris coworkers 
have employed vortex methods for the study of turbulent combusting flows [131,132]. Beljaars et al. 

[133] have developed a structural model for turbulent wall bounded flows in which they solve the 
viscous flow equations for the inner flow and solve the outer flow using point vortex method. Dai et al. 

[134] have studied the flow over a backward facing step using Chorin’s random vortex method. They 
found that the vortex method predicted the essential characteristics of internal recirculating turbulent 
flow. The calculated reattachment length agreed well with the experimental results. The equations of 
motion of point vortices represent a Hamiltonian system. The study of the dynamics of few vortex 
systems is expected to shed some light on the transition from laminar to turbulent flow and on the onset of 
chaos [126]. 


9.2 Large-Eddy Simulation 

In large-eddy simulation, the filtered time dependent three-dimensional Navier-Stokes equation is 
integrated directly. Filtering is done to remove the small eddies and to obtain an equation for large 
eddies. It is essentially equal to averaging over a small spatial region or low-pass filtering the equations 
in Fourier-space. The effect of small eddies on the large ones has to be modeled. For a given flow, the 
dividing line between the small and large scale flow fields is determined by the computational grid reso- 
lution. When a large number of grid points is used, a large fraction of the turbulent eddies is direetly 
calculated and only a small fraction has to be modeled. One important feature of LES is that most of 
the time-dependent characteristics of the flow are retained rather than lost in the averaging process. 
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To obtain the equations for LES, Leonard [135] suggested the local spatial averaging: 


Uj(x) = J" G (x-x') u- (x') dx' (58) 

where G is a normalized weighting function, which may be called a filter. Its effect is to remove the small 
scale fluctuations from Uj in forming u“ the large scale field, and Uj', the subgrid scale field. A simple 

choice is to let G = 1 within a cubic volume with sides of length centered at x and let G = 0 outside. 
Then, 


Uj(x) = 



x^-x^ , x^-x^') dxj' dx2* dx-j' 


(59) 


The integrations being from Xj-1/2 A a to Xj + 1/2 A a , x-> - 1/2 A & to X 2 + 1/2 A a and X 3 - 1/2 A a 
to x^ + 1/2 A g . 

The filtered counterparts of (1) and (2) can be obtained by multiplying (1) and (2) by the 
weighting function G = 1 and integrating over the cubic volume v 


dUj 



(60) 


3u- 

"at" 


+ 




dp 

3xj 


+ v V 



(61) 


Now make the substitution u- = u- + u-' in the nonlinear advective term to obtain 


!!!l + JL 

aT + 



3p 9 __ 9 

_ + „V Uj- — % 


(62) 


where 


Tjy = U iUj 


+ Uj Uj + U 


i 11 ] 


(63) 


is the subgrid scale (SGS) Reynolds stress. Thus in LES one solves equations (60) and (62) together with 
(63). Since the small scale component of the velocity field Uj' is not computed, the terms containing it 
need to be modeled. 


52 



Deardorff [136] and Schumann [137] employ a different approach to obtain LES equations. 
The equations are set up for the grid employed based on integral conservation equations for each grid 
volume. The form of the equations are the same as (60) and (62) if the overbar is interpreted as the 
volume average. In this approach 


u^Uj = Uj Uj . (64) 

^ = 0 • (65) 


as in the time-averaging process. This reduces the SGS Reynolds stress to 


T ?ij = u 



( 66 ) 


Uj'uj' is modeled by eddy viscosity models. 

The model often used is that due to Smagorinsky [138] , 



where C is a consta nt wh ich varies from 0.13 to 0.21. [When using the form (63), the term to be 
modeled is still only u^uj as the other terms can be taken care of by using equation (64) in equation 

(61) which results in a modified pressure gradient term [139] ]. 

The near wall treatment is crucial as in the time averaged Reynolds stress modeling. For a success- 
ful LES simulation, one has to be able to resolve all the large eddies in the neighborhood of the wall 

and at a distance from the wall, since the near wall eddies are responsible for the turbulent energy pro- 
duction for the flow field. Moin and Kim [140] solve the LES equations up to the wall while Schumann 
[137] and Deardorff [136] carry out the calculations only to a point in the logarithmic layer where 
semi-empirical boundary conditions are used to near wall turbulence. In this approach the important 
near wall turbulence flow dynamics is modeled which limits the use of LES approach to the understand- 
ing of the flow physics. 

The channel flow has been computed using the LES approach and found to reproduce experi- 
mental observations such as the turbulent bursting sequence, production of longitudinal vortices in 
boundary layers [138-140]. An example of the LES solution is shown in Figure 34. The figure shows 
contours of instantaneous velocity, pressure fluctuations and SGS kinetic energy in a channel flow. Clark 
et al. [141] evaluated the SGS Reynolds stress model constant and found that it agreed with the one 
obtained from theory and experiment. Though it is expected that the LES and direct simulations can help to 
improve the time averaged Reynolds stress modeling [139], no such improvement has ensued. 
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Figure 4. Contours of velocity, SGS kinetic energy and pressure fluctuations in a channel flow [137], 






9.3 Direct or Full Simulation 


In this approach one solves the time dependent three-dimensional Navier-Stokes equations for 
different scales of the flow (the smallest scale that can be resolved being determined by the number of 
grid points). This approach has been used mainly to homogeneous flows at low Reynolds numbers (of 
the order of R t = 40), employing as much as 2 million grid points. Periodicity is assumed in all three 

dimensions. The velocity field is represented at every instant by a three-dimensional Fourier series in a 
coordinate system moving with the spatially linear mean flow [142]. In these coordinates the full 
Navier-Stokes equations admit spatially periodic solutions with fixed period for all time. The representa- 
tion of statistical homogeneity in space by strict periodicity requires that the computational period be 
much longer than any turbulence scale containing significant energy. Also, the mesh size has to be much 
smaller than the turbulence dissipation scales. A low turbulence Reynolds number is required to limit 
the scales. With the present computer resources, only flows with a small range of energetic turbulence 
scales can be computed [142]. 

Large eddy simulation and direct simulation are not expected to become engineering calculation 
methods in the foreseeable future. However, these methods will help in the understanding of the physics 
of turbulent flows. 


X. CONCLUSIONS 


The report presents a brief account of various turbulence closure models and their applications 
to internal flows. A critical evaluation of the performance of the models is presented. The main con- 
clusions are: 

1 . The k-e model is the most widely used model for internal flows. It provides an efficient 
method of calculating engineering flows. The performance of the standard k-e model becomes poor as we 
go from the attached flow, recirculating flow, swirl flow, to combusting flows in that order. 

2. The standard k-e model performs poorly in the prediction of apparently simple flow over a 
backward facing step. There is no significant improvement in the value of reattachment length predicted 
with modified near-wall models. Algebraic stress model should be used with modified e-equation. 


3. The k-e model performance is rather poor for flows with streamline curvature. Algebraic stress 
model performs better for these flows. 

4. Reynolds stress models have not been thoroughly tested for recirculating flows and swirl 
flows. Computational effort required for RSM is much greater than that required for the k-e model. 

5. For flows with regions of secondary flow (noncircular duct flows) algebraic stress model 
performs fairly well for fully developed flow. For developing flow RSM should be employed. 

6. In many computations, the grid dependence of the solutions has not been achieved or tested 
due to the limitations of the computing resources. These solutions have to be taken with caution. 
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7. The initial conditions play a crucial role in the performance of the model and predicted 

solutions. Thus, use of appropriate initial conditions in the SSME flow computation is essential. 

8. The TEACH (Teaching Elliptic Axisymmetric Characteristics Heuristically) code [143] or its 
variant is the most widely used computer code in the prediction of turbulent internal (plane 2D, 
axisymmetric, 3-D, swirl and reacting) flows, cited in the literature. 

XI. RECOMMENDATIONS 

Based on the present evaluation, the following suggestions for future work can be made: 

1. More work is necessary to arrive at a better numerical scheme for the solution of the model 

equations. 

2. A clear statement of the initial conditions employed should be mandatory for any presenta- 
tion of the model predictions, especially when they are made “suitable” to match the experimental 
results. 


3. The use of ASM models for flows with streamline curvature appears promising and needs 
further testing. 

4. The use of the flux Richardson number for curved flows should be explored. 

5. Detailed experimental data for curved flows with strong curvature are needed for the valida- 
tion of the model predictions for these flows. 

6. With the availability of more computing power, use of RSM models to account for curvature 
effects, countergradient transport and secondary flows would improve the confidence in turbulence 
closure models. 
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